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Charging of a clean two-dimensional island is studied in the regime of small concentration of 
electrons when they form the Wigner crystal. Two forms of electron-electron interaction potential 
are studied: the pure Coulomb interaction and the exponential interaction corresponding to the 
screening by a pair of close metallic gates. The electrons are assumed to reside in a parabolic 
external confining potential. Due to the crystalline symmetry the center of the confinement can be 
situated at distinct positions with respect to the crystal. With the increasing number of electrons A*" 
the center periodically hops from one such a location to another providing the lowest total energy. 
These events occur with the period ~ N^^"^ . At these moments in the case of the pure Coulomb 
interaction the charging energy of the island has a negative correction « 15%. For the case of the 
exponential interaction at the moments of switching the capacitance becomes negative and ~ A'^^'^'* 
new electrons enter the island simultaneously. The configurations of disclinations and dislocations 
in the island are also studied. 

PACS numbers; 73.20.Dx, 73.40. Gk, 73.40. Sx 



(N 
> 

O 
(N 
l> 
O 



I 

o 
o 



I. INTRODUCTION 

In recent experinientsl'i the charging of a quantum 
dot is studied by the single electron capacitance spec- 
troscopy method. The quantum dot is located between 
two capacitor plates: a metallic gate and a heavily doped 
GaAs layer. Tunneling between the dot and the heav- 
ily doped side is possible during the experimental times 
while the barrier to the metal is completely insulating. 
DC potential Vg and a weak AC potential are applied to 
the capacitor. With the increase of Vg the differential 
capacitance experiences periodic peaks when addition of 
a new electron to the dot becomes possible. The spac- 
ing between two nearest peaks AVg can be related to the 
ground state energy E{N) of the dot with A^ electrons: 



aeAVg = E{N + 1) - 2E{N) + E{N - 1) 



A(A^) = eye 



(1.1) 



Here a is a geometrical coefficient, A(A^) is the charging 
energy, Cn is the capacitance of the dot with A^ elec- 
trons. It was observed in Refs. |^,^ that at a low con- 
centration of electrons or in a strong magnetic field the 
nearest peaks can merge, indicating that at some values 
of Vg two or even three electrons enter the dot simultane- 
ously. In other words some charging energies apparently 
become zero or negative. In a fixed magnetic field this 
puzzling event repeats periodically in A^. Disappearance 
of the charging energy looks like a result of an unknown 
attraction between electrons and represents a real chal- 
lenge for theory. 

Pairing of the differential capacitance peaks has been 
studied theoretically before for disordered dots. Explana- 
tion of the pairing based on the lattice polaronic mech- 
anism has been suggested in Ref. ||. In Ref. ^ it was 



demonstrated how electron-electron repulsion, screened 
by a close metallic gate, can lead to electron pairing for 
a specially arranged compact clusters of localized states 
in a disordered dot. This effect is a result of redistribu- 
tion of the other electrons after arrival of new ones. It 
was interpreted in Ref. ^ as electronic bipolaron. 

In this paper we study the addition spectrum of a dot 
in which the density of electrons is small and the exter- 
nal disorder potential is very weak, so that electrons in 
the dot form the Wigner crystal. We call such a dot a 
Wigner crystal island. In the experimental conditions of 
Ref. I one can think about a Wigner crystal island lit- 
erally only in the highest magnetic field. One can also 
imagine similar experiments with a Wigner crystal island 
on the surface of liquid helium. In the present paper we 
consider the extreme classical limit of the Wigner crys- 
tal, when the amplitude of the quantum fluctuations is 
much smaller than the interparticle distance. In this case 
one can think of electrons as of classical particles and the 
energy of the system is given by the following expression: 



(1.2) 



i<3 



The first term represents interactions among electrons 
located at points r^, with U {r) being the interaction 
potential. The second term is the contribution to the en- 
ergy due to the external confinement, which is assumed 
to have parabolic form. Coefficient A plays the role of 
strength of the confinement. The forms of the interaction 
potential considered are: the pure Coulomb interaction 



U (r) — e^ / K.r, 



(1.3) 



with e and k being the electron charge and the dielectric 
constant correspondingly, and the exponential interac- 
tion 
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U (r) = t/o exp {—r/s) . 



(1.4) 



The latter interaction potential corresponds to the case 
when the island is situated between two metallic gates, 
with s being of the order of the distance between them 
and Uq ^ e^/ns. In this paper we study the case a:s> s, 
where a is the average interparticle distance. 

First we study the addition spectrum E{N) of such a 
system numerically. We notice that for the both types of 
interaction the energy of the ground state has a quasi- 
periodic correction. The period of this correction scales 
as ~ y/N, or the number of crystalline rows in the is- 
land. Its shape is universal, i.e. independent of the form 
of the electron-electron interaction. Explanation of these 
oscillations is the subject of the subsequent theoretical 
analysis. 

We attribute the aforementioned oscillations to the 
combination of two effects: insertion of new crystalline 
rows and motion of the center of the confinement rela- 
tive to the crystal. Consider the former effect first. Let 
us fix the position of the center of the external parabola 
relative to the adjacent Bravais unit cell. For example let 
it coincide with one of the lattice sites. As electrons are 
added to such a system, the number of crystalline rows 
grows roughly as ^ ^/N. The periodic appearances of 
the new rows bring about oscillations of the total energy 
E{N) with N. The period of these oscillations scales as 
SN ^ \/N in agreement with our numerical results. 

Let us discuss the influence of the position of the con- 
finement center. The energy of the system as a function 
of this position has the same symmetry as the lattice. 
Hence the extrema of the energy have to be situated at 
the points of high symmetry, e.g. the centers of two- 
fold rotational symmetry or higher. There are three such 
points in the two-dimensional (2D) triangular lattice (see 
Fig. ||a below) . The evolution of the island with N con- 
sists in switching between these three positions of the 
center, every time choosing the location providing the 
lowest total energy. This effect is only based on the sym- 
metry considerations and is therefore universal, i.e. in- 
dependent of the form of interaction. This idea suggests 
a simple recipe for the calculation of the energy of the 
island. One has to calculate three energy branches corre- 
sponding to the different locations of the center and then 
choose the lowest one. 

Let us consider these two effects in more detail, sep- 
arately for the two types of interaction. We start with 
the case of extremely short-range interaction given by 
Eq. (1.4): a ^ s. This case can be realized, e. g., if the 
confining parabola is very weak compared to the interac- 
tion prefactor: AR^ <C Uq. Here R is the radius of the 
island. For this case the interaction between the nearest 
neighbors is a very steep function of distance due to the 
fast decay of the exponential (L4). The variations of the 
lattice constant in this case are very small (see a more 
elaborate discussion in Section [II). This implies that 
the crystalline rows are almost straight lines (see Fig. |^) . 
Hence we deal with a piece of almost perfect triangular 



crystal with new electrons being added on the surface of 
the island. A new crystalline row therefore appears on 
the surface. It can be shown that this creates an anoma- 
lous increase in the density of states (DOS) of electrons. 
Such variations of DOS result in the appearance of the 
periodic correction to the energy of the island. This cor- 
rection is considered in detail in Section |l| An interest- 
ing implication of this picture is the multiple electron en- 
tering. It turns out that if one slowly raises the chemical 
potential of electrons in the island, then at the points of 
switching between the branches mentioned above about 
/\ri/4 electrons enter the island simultaneously. A simple 
model of this phenomenon has been suggested before in 
Ref. H for a small island containing < 55 electrons. In 
this paper we discuss this phenomenon for larger islands. 
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FIG. 1. The configuration of electrons in the island 
for the case N^80,A = 10~®, (7o = 1, and s = 1. A typ- 
ical surface particle in a large dot has four nearest neigh- 
bors. Two particles with three and five nearest neighbors 
are shown by the circle and the triangle correspondingly. 
They represent a pair of opposite surface disclinations or 
a surface dislocation. 

One can think about the new rows appearing in the 
crystalline island as of pairs of dislocations of opposite 
sign. In the case of the short-range interaction these 
defects are pushed to the surface by a huge price for 
elastic defo rmations. In the Coulomb case, described 
by Eq. (|l.3|), the shear modulus and, subsequently, the 
Young's modulus of the crystal are relatively small. As 
discussed in Section ^ the incommensurability of the cir- 
cular shape of the dot with the lattice and the inhomo- 
geneity of the density of electrons generate in this case 
topological defects, disclinations and dislocations, inside 
the island. We argue that these defects determine the 
variations of the energy when the center of the confine- 
ment is fixed relative to the crystal. The number of dis- 
locations scales as the number of crystalline rows ^ \/N. 
This implies that a new dislocation appears every ^ ^/N 
electrons. Due to the discreteness of these defects a single 
branch of the total energy corresponding to a fixed po- 
sition of the center acquires a quasi-periodic correction. 
This periodic building up and relaxation of the elastic 
energy of the dot caused by the discreteness of disloca- 



2 



tions is very similar to the variations of the electrostatic 
energy brought about by the discreteness of electrons, 
known as Coulomb blockade. Following this analogy we 
call the former periodic phenomenon an elastic blockade. 

The elastic blockade appears to be a bit more compli- 
cated than its electrostatic counterpart. The center of 
the confinement can move among three distinct points of 
the triangular lattice mentioned above making the sys- 
tem switch from one energy branch to another. Such a 
sudden switching results in a correction to the charging 
energy [see Eq. (1.1)]: 



SA 



-0.15A, 



(1.5) 



where A is the average charging energy. This reduction 
of the distance between the nearest Coulomb blockade 
peaks happens with the period ~ -s/iV determined by the 
elastic blockade. It can be thought of as an analog of 
merging of a few peaks observed in the case of the short- 
range interaction. 

The optimum configurations of electrons in the 
parabolic confinement with the Coulomb interaction have 
been studied before in Ref. ^. Our numerical results both 
for the energies and the configurations agree with the re- 
sults of this work. However our interpretation of the re- 
sults is different. The authors of Ref. H adopt a model in 
which the electrons fill in shells concentric to the perime- 
ter of the island. We argue that in the regime studied 
in the numerical experiment only a narrow ring adjacent 
to the perimeter is concentric to it. The width of such a 
ring is ~ y/ R ■ a. The rest of the island is filled with an 
almost perfect crystal, (see a Sections ^ and [V^ . 

The paper is organized as follows. First we consid er th e 
case of an extremely short-range interaction [Eq. ( |1.4D ], 
when the Young's modulus of the crystal is very large and 
no lattice defects can exist in the interior of the island. 
This case is discussed in Sections || and In Sec- 
tion IV we report the results of the numerical sol utio n 



of the problem with the Coulomb interaction [Eq. (1.3)]. 
In Section ^ we turn to the discussion of different kinds 
of defects, which can exist in a compressible island. In 
Section we discuss the theory of the elastic blockade 
for the case of the Coulomb interaction. Section VI] is 
dedicated to our conclusions. 



II. SHORT-RANGE INTERACTION: 
NUMERICAL RESULTS 

In this section we consider the system of N electrons 
interacting by a strongly screened short-range potential. 
This limit is defined by Eq. (1.2) with interaction given by 
Eq. (1.4) and a ^ s. As it is shown later in Section III the 
latter condition can be realized if the interaction prefac- 
tor significantly exceeds the typical confinement energy: 
Uq ^ AR^. At this condition the variations of the lattice 
constant of the crystal are negligible. This indeed can be 
observed in the configurations obtained in the numerical 



experiment. One of such a configurations is shown in 
Fig.|l|, obtained at a « 13s. 

In our numerical analysis we minimized the energy 
fmictional (1.2) with the exponential interaction given by 
Eq. (1.4) with respect to the coordinates of electrons us- 
ing a genetic algorithm similar to that outlined in Ref. |^. 
Below we describe this numerical technique and the mo- 
tivation for its use. 

The problem at hand belongs to the vast class of 
problems of finding the global minimum of a multi- 
dimensional function, which has plenty of local minima. 
The well developed methods for convex differentiable 
functions (the conjugate gradient, Newton's method) do 
not work here as they find only some local minimum. The 
most frequently used method in thi|S| case is the Metropo- 
lis simulated annealing technique.El In this method the 
system is modeled at some artificially introduced tem- 
perature, which is gradually decreased to a very small 
value. It is assumed that after this annealing the sys- 
tem falls into the state of the lowest energy. Although 
in principle in this method the system can hop from a 
metastable state to the ground state, in practice, if the 
potential barrier between them is high, the time to per- 
form such a hop can exceed the time of the simulation. 

An alternative method is the genetic algorithm, which, 
was proved to be superior to the simulated annealing.LJ 
Initially five different parent configurations were obtained 
by relaxing random configurations using the conjugate 
gradient algorithm. Then these configuration were mated 
pairwise to obtain additional 15 child configurations (in- 
cluding mating with itself) , which were again relaxed us- 
ing the conjugate gradient algorithm. Mating consisted 
in cutting two configurations into halves by a random 
line and then connecting those halves belonging to dif- 
ferent parents to form a new child. From the resulting 20 
conformations (parents and children) five were chosen to 
be parents for the next iteration. The new parents con- 
sisted of the lowest energy configuration and the other 
four separated by at least the precision of the conjugate 
gradient. This has been done to avoid dominating the 
process by one conformation. During ten iterations as 
many as 155 local minima were examined. The energies 
of the optimum conformations obtained in this way agree 
with those of Ref. ^ and for some N are even lower. Be- 
fore presenting these results we would like to discuss the 
method we used to process the data. 

To analyze the dependence of the ground state energy 
E (N) on the number of particles we split it into the 
smooth E (N) and the fluctuating component 5E (N) : 



E{N)^ E (N) + SE {N) 



(2.1) 



in the manner of Ref. ^. The smooth component has the 
form E{N) = r]iN^ + 772 + mN"^^^ + r]4N^/^ + r]^N^/'^, 
where the coefficients f?i • • • ?75 are chosen to minimize the 
fluctuations. The fluctuating part is displayed in Fig. ||. 
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FIG. 2. The fluctuating part of the ground state en- 
ergy in the model with the exponential interaction ( |l.4| ) 
for A = 10~^ Uo = 1, and s = 1. 

The curve in Fig. |^ evidently has a quasi-periodic 
structure. It consists of the series of interchanging deep 
and shallow minima, separated by the peeks with more or 
less smooth slopes. We call this structure the "Kremlin 
wall" for its resemblance to the embattlements on top of 
the walls of the Moscow Kremlin. Its period with a very 
good precision scales as N^/"^ and numerically is equal to 
the number of electrons in the outer crystalline row of 
the island. The amplitude of the oscillations grows with 
N afi 



\6E\ozN'', 7 = 0.8 ±0.1. 



(2.2) 



Although the crystal in the bulk is almost perfect, its 
boundary is extremely irregular (see Fig. |l|). It can be 
thought of as a superposition of various types of defects 
pushed against the surface by extremely large Young's 
modulus, provided by the short-range interaction. These 
defects can be associated with the particles having an 
anomalous coordination number. Normally a particle 
on the surface of the triangular crystal has four near- 
est neighbors. But there are particles having a coor- 
dination number equal to three or five. These particles 
can be associated with positive and negative disclinations 
on the surface correspondingly, as their creation assumes 
removal (insertion) of a 7r/3 wedge (see Fig. ||). Dislo- 
cations on the surface are pairs of such disclinations of 
the opposite sign forming a dipole. At a small Young's 
modulus however these defects can dive inside the island. 
This transition is quantitatively described in Section 0. 



III. SHORT-RANGE INTERACTION: A THEORY 

In this section we use a hard disk model to explain our 
numerical results obtained in the limit Uq ^ AB? . To 
justify this model we first show that the variation of the 
lattice constant in the island is small in this limit: 

,5a = a(i?) - a(0) = ^ In TV < a(0). (3.1) 

Then we demonstrate that the interaction energy is small 
compared to the confinement energy and hence can be 
neglected. 



To prove the first statement we find the pressure a„ 
in the crystal associated with its contraction by the ex- 
ternal potential Ar^ . The solution can be found similar 
to Ref. i 



CT,, = -S[(t)A[R 



2 2\ 

r 



(3.2) 



where S{a) = (3 + cr)/4, 3/4 < S{a) < 1, and ct is the 
Poisson ratio. The solution is easy to understand for the 
limiting case of liquid cr ^ 1 , when Oik = —pSik , p being 
pressure. In this case p is of purely hydrostatic origin: 



dp{r) 
dr 



-2Arn{r). 



(3.3) 



Assuming the density n(r) to be uniform and p{R) = 
we obtain the following solution: 



p{r) = AniR"^ - r^), 



(3.4) 



which agrees with Eq. (3^) taken at cr = 1. The forces 
produced by this pressure have to be balanced by the 
interaction forces between particles: 



An{R 



2 2\ 

r a ' 



(3.5) 



In this equation we assumed S{(7) ^ 1. From this condi- 
tion we obtain the following result for the lattice spacing 
in the island: 



a(r) 



;ln 



Un 



A{R^-r^) 



(3.6) 



The difference between the lattice spacing on the bound- 
ary and in the center can be readily obtained: 



R 



(5a w s In ( — — ^ I — s In I —At 1 = s In 
ARaJ \AR^ I a 



In N. 



(3.7) 



Comparing ( |3.6| ) and (3/7) one can see that the case 
a ^ 6a can be realized when Uq ^ AR^. This im- 
plies that the prefactor of the interaction potential has 
to be large. In this case the nearest particles are far 
away from each other where the exponential interaction 
is very steep. This eliminates significant changes in the 
interparticle distance. 

Although the derivative of the interaction potential is 
important in determining the lattice spacing the typical 
value of interaction energy in this case is very small. The 
characteristic interaction energy per particle ca n be es- 
timated as: Enn ~ Uoexp{—a/s). From Eq. ( |3.5| ) we 
conclude that 



a 



(3.8) 



Hence in the considered regime a 3> s the interaction 
energy indeed can be neglected. 
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It is plausible then to accept the hard disk model to 
explain our numerical results obtained for this case. We 
assume therefore that the interaction has the form: 



Uir) 




(3.9) 



In this model the interaction energy is zero and the total 
energy of the system can be written as: 



N 
i=l 



(3.10) 



where belong to the triangular lattice, with the lattice 
spacing equal to the radius of the interaction a. This en- 
ergy formally coincides with the moment of inertia of the 
system of N particles of mass A. Hence to find the mini- 
mum energy configuration one has to cut a piece from the 
triangular crystal that has a minimum moment of inertia. 
This piece must contain the given number of particles A''. 
The average values of the total energy, chemical potential 
and the charging energy are given by: 



E^AN- = AN^/2Trn 
p, = dE/dN = AN/nn 
A = dfl/dN = A/nn, 



(3.11) 



where R is the average radius of the circle, mrR^ 
n = 2/a^V3 is the concentration of lattice sites. 

A similar problem has been considered in the litera- 
ture in the context of the so-called circle problem. The 
problem is to find the fluctuations of the number of par- 
ticles belonging to some lattice inside of a circle of radius 
R (see Ref. ^). The quantity which has been studied is 
the deviation of the number of particles from its average 
value: 



SN{R) = N{R) - nnR^ 



(3.12) 



where mrR^ = N is the average number of particles in 
the circle. The classical circle problem which goes back 
to Gauss is to find the uniforrn bound for \SN{R)\. The 
best result in this direction iaj 



\SN\ < C,R^^^'^+' ^ iv23/73+^/2 



(3.13) 



Here is a i?-independent constant and e > 0. A sim- 
ple estimate for these fluctuations can be obtained from 
the assuroption that they occur around the perimeter of 
the circleH in the strip of width of the order of the lattice 
spacing a. The number of such points A'pcr ^ R/a and 
the fluctuations are 



1/4 



(3.14) 



This estimate may imply a pseudo-random behavior of 
SN. N otice that it agrees with the uniform upper bound 
(^l3|) as 23/73 > 1/4. 



Ref. ^ also considers the distribution of SN. It was 
found that the distribution function is non-Gaussian: 



p{SN) < exp {-6N^/N) , 



(3.15) 



with the mean square deviation given by Eq. (3.14). It 
was also noted that 6N{R) in addition to be pseudo- 
random is almost periodic function of the radius. 





I ~ (Ra)'" 
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FIG. 3. The circular island (gray) covering the crys- 
talline rows, shown by lines. The lattice points within the 
circle are shown by larger disks than those outside. Every 
time the circle intersects a new row a new terrace appears. 
One of such terraces is shown by the bold segment. 

Below we will concentrate on the periodicity of 6N{R), 
rather than on its randomness. In essence we would like 
to find the average of dN{R) over many periods and un- 
derstand the nature of its oscillations. Then the question 
about the deviations from this average becomes relevant. 
Our analysis shows that the oscillations are associated 
with the periodic intersections of the circle with the crys- 
talline rows (see Fig. ^). At the moment of intersection 
an anomalously large number of particles can be put into 
the circle. This brings about an increase in the func- 
tion SN{R) every time the circle hits a new crystalline 
row. This increase is accompanied by a decrease later in 
the period as the average value of SN is zero. To illus- 
trate this point consider the triangular lattice (Fig. |^). 
In this case for the most part the oscillations of SN are 
associated with the collisions of the circle with the main 
crystalline rows, {VS, 1) plus other five obtained by 7r/3 
rotations, separated hy h = Hence the variations 

of SN associated with these rows have a period in R equal 
to Tfi = h. One can collapse SN from many such periods 
into one and average over all of them. To do this we first 
define a normalized variation of the number of particles: 



(3.16) 



According to Eq. (3.14) the average amplitude of this 
function does not change with R. To extract the com- 
ponent of this function having the period Tr — h we 
average it over many such periods: 



Srj (R) = lim 



1 

— ^ 577 (TRm + R) 



m=0 



(3.17) 
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Thus obtained function SN{R) = Srj (R) is shown in 
Fig. ^ by a smooth hne. The actual 6N{R) is also pre- 
sented in this figure to show that its main variation is 
indeed associated with the aforementioned period. 




FIG. 4. The averaged oscillations of the number of 
particles as a function of the radius of the circle (bold 
smooth line) and the actual value, observed in the numer- 
ical simulation on the lattice (thin rough line). The center 
of the circle is located at one of the lattice sites. 

Before embarking into the calculation of the 6N we 
would like to mention that other crystalline planes pro- 
duce similar oscillations. These oscillations have a period 
that is smaller than h and can be incommensurate to it. 
The amplitude of these satellite oscillations is numeri- 
cally smaller than the main one and hence they will be 
neglected in this paper. All the periodicities added to- 
gether produce a complex fractal curve shown in Fig. ^. 

Before getting into details of the averaging we would 
like to show a simple way to estimate the amplitude of 
these oscillations. Look at Fig. ^. The collision of a crys- 
talline row with the circle results in the formation of the 
"terrace" on the surface of the crystal bounded by the 
circle. The average length of such a "terrace" is 



I - VR-a. 



(3.18) 



The number of crystalline sites in such a "terrace" pro- 
vides an estimate for the variations in the number of 
particles inside the circle: 



SN = l/a - y/R/c 



(3.19) 



This estimate agrees with Eq. ( ^.14 ). However, it is 
based on the existence of the ordered periodic structure 
in SN. 

This argument allows us to estima te the amplitude of 
the oscillations of the total energy ( |3l0| ). Indeed, we 
have found the oscillations of the number of particles as 
a function of radius or, in other words, of the chemical 
potential fi = AR^. Now it easy to calculate the related 
oscillation of the chemical potential as a function of the 
number of particles. To this end we notice that if the os- 
cillations are weak they are simply proportional to each 
other: 



Sn = -SN/D, 



(3.20) 



where = 1/A is the average density of states. The pe- 
riod of S^{N) is related to the period of SN{R) in the 
obvious way: Tn = dN/dR ■ Tr = nnRh ^ Now 
the oscillating part of the total energy can be easily esti- 
mated: 



SE{N) = / dNS^i{N) 

- Tn ■ SNA - Aa'^N^/^. 



(3.21) 



Notice that both the periodicity and amplitude of the os- 
cillation s gi ven by this estimate agree with our numerical 
results (2.2). Below we will calculate the form of these 
oscillations averaged over many periods. 

We now turn to the calculation of SN{R). We follow 
the guidelines of Ref. ^. Assume that the center of the 
circle is located at point Tq reduced to Bravais unit cell. 
The number of particles in the circle can be expressed 
through the sum over the lattice points Ri . 



6N{R) ^ N - 7171 R^ 



Y^fiRi)- [ nd\f{r), 
Ri 



(3.22) 



where 



1, |r-ro| <i?, 
0, |r-ro| > R. 



(3.23) 



Using the Poisson summation formula, Eq. (3.22) can be 
rewritten as follows: 



SN[R) 



(3.24) 



where the summation is assumed over the reciprocal lat- 
tice vectors and / (q) = 2-KRJx{\ q\ K) l |<7|-exp(— iqro) 
is the Fourier transform of function (B.23), with Ji{x) be- 
ing the Bessel function. For the normalized variation of 



number of particles (B.16) we then obtain 



(3.25) 



The hmit in Eq. (|]T|) can be most easily evaluated 
upon noticing that Ji(x) « ^/2/^Tx cos {x — 37r/4), when 
X ^ 1. Changing the order of summation over Qi and 
m we obtain: 



Sr] (R) = lim 



M- 



(3.26) 



^ cos (I I Tnm +\Q,\R- 37r/4) ^_,Q^,„ 
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To evaluate the sum over m we again can use the Poisson 
summation formula, which establishes a selection rule for 
the values of the wave- vector Qf. 



\Q,\TR^27rn, 



(3.27) 



where n is an integer. The reciprocal lattice vectors are 
representable in the form Qi = Qq (fci-\/3/2, fci/2 + fc2), 
where Qq = Air/ay/S and ki^2 are integers. Equation 



( |3.27 ) can be rewritten in term of ki 2 in the following 



way: 3fc^ + (fci + 2/02)^ = 4n^, which is a Diophantine 
equation. It has a trivial set of solutions, ki — 0, k2 = n 
plus five others obtained from it by 7r/3 rotations. They 
correspond to the collisions of the circle with the main 
sequence of the crystalline rows, separated by Tr. These 
solutions will be taken into account below. The other so- 



lutions of Eq. ( 3.27 ) are more sparse in the Q-space, and 
hence produce smaller than Tr periods in the real space. 
They correspond to the collisions with other crystalline 
lines. They produce numerically smaller variations than 
the main sequence and will be disregarded below. Thus 
we conclude that 



lim 

M^oc 



1 

— J2 cos imTRTn +mR- 3TT/4) 



m=0 



Y.^Q..Q. cos(|Q,|i?-3V4), 



(3.28) 



where Qi are solutions of (3.27) with the reservations 
mentioned above. The remaining summation over Qi is 
easily performed and we obtain finally the result for the 
average variation of the number of particles: 



(3.29) 



where are the unitary vectors normal to the six main 
crystalline rows, e,; = t/* (a/3/2, 1/2), IJ is the 7r/3 ro- 
tation, i = 0, ... ,5. Function 6Nq determines the oscil- 
lations produced by one set of rows if the center of the 
circle Tq coincides with one of the lattice sites: 



23/433/8 

5No (R) = 



ri/4 



1^1- (»») 



Here 



2r(l - z) ^ sin (27rfcg -h z7r/2) 



(3.31) 



fc=i 



is the generalized Riemann zeta-f uncti onB In addition 
to the numerical coefficient in Eq. ( 3.3Cl| ) we obtained the 
overall amplitude N^^'^, which agrees with our qualitative 
analysis, and the oscillating part {—1/2, R/Tji), which 
has an amplitude of the order of unity. In Fig. ^ 5 Nrn=o 
is shown by a smooth line. It is clear from Eq. ( 3.29| ) 



that six sequences of the crystalline rows contribute to 
oscillations independently, with e^ro determining a "de- 
lay" produced by the displacement of the center of the 
circle from a lattice site. 

We would like now to calculate the oscillations of the 
total energy of the incompressible island. First, as pa- 
rameters of the problem we will use Tq and R. This 
implies that we will have a given value of the chemical 
potential in the island /i = AR^ . The relevant thermody- 
namic potential in this case is fir oi^J■)■ We will evaluate 
this potential and optimize it with respect to tq. We will 
then use the theorem of small increments to find E{N). 
This function will indeed be similar to the Kremlin wall 
structure displayed in Fig. ^. 

The first step in this program is realized similar to ob- 
taining 5N. We apply the method used above to the 
function 

n{R) = E{R) - fiN, n = AR^. (3.32) 
The average value of this function is negative: 

n = ~AN-R^/2^~E, (3.33) 



(compare to Eq. (3.11)). The deviation of the potential 
from the average is given by 

6il{R) ^6E{R)~ AR^6N{R) 

J29{Ri)- [ nd\g{r)~AR^5N{R), 

(3.34) 



R. 



where 



7(r) 



A(r-ro)^ |r-ro| < R, 



0, 



|r-ro| > R. 



(3.35) 



Using 



Ji(|g| R) 2 — J2(|g| R) 



(3.36) 



one can obtain the oscillations of the thermodynamic po- 
tential in exactly the same manner as those of the num- 
ber of particles. To this end, as it follows from ( 3.21 ), 
one has to average the normalize d va riations of poten- 
tial: 5uj = SVt/R^/"^ (compare to ( |3.16| )). Eventually we 
obtain 



5nr„ {R) = '^^0 - ^'^0) 



(3.37) 



where 



m-,{R) = -A'-^^[N{R)f\[-l,^^, (3.38) 
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is the function describing the oscillations if the center of 
the circle coincides with one of the la ttice sites. Its am- 
plitude agrees with o ur presumption (3.21). Notice that 
Sflro given by ( 3.37 ) is in fact a function of the chem- 
ical potential as the latter is related to the radius by 
/i = AR^. This function for rp = is shown in Fig. ||. 




FIG. 5. The averaged oscillations of the thermody- 
namic potential f2 as a function of the radius of the circle 
(bold smooth line) and the values actually observed in the 
numerical simulation (thin rough line) . The center of the 
circle is located at one of the lattice sites. 

The next step is to minimize it with respect to Tq. To 
do this we notice that Sflrg as a function of ro must have 
the same symmetry as the lattice. Hence its extrema 
have to be located at points denoted in Fig. ||a by A, B, 
and C, unaffected by the transformations of the symme- 
try. Those are the points coinciding with a lattice site, 
the center of the triangular face, and the point half-way 
between two lattice sites. Thus we actually need to make 
a choice between these three positions of the center. This 
choice has to be made to minimize the energy SQro- In 
Fig. ^3 this minimization is done graphically. From the 
three curves S^Ia, S^Ib, and 6flc for every value of R the 
lowest one is chosen. The resulting curve is shown by a 
bold line. It indeed looks like the Kremlin wall structure 
obtained in the numerical experiment. 

We would like now to discuss the points of switching 
between two curves. At these points the center of the 
island hops from one of the locations A, B, or C to an- 
other. If the chemical potential is below or above one of 
these points the newly added electrons are spread uni- 
formly over the circumference of the island, so that the 
center of mass does not move. But exactly at this point 
~ N^^^ electrons travel from one side of the island to an- 
other shifting the center of mass relative to the crystal. 
In addition to that at the point of transition ^ N^^^ new 
electrons enter the island. To understand this we recall 
that TV = —dil/d^. As it is seen from Fig. ^ at the point 
of transition the first derivative of has a discontinuity 
of the order of 



A 



dsn 

dfi 



1 



2AR 



-A 



dsn 

'dR 



_Aa^N-^ 



(3.39) 



This implies that at this point a few peaks of the 
Coulomb blockade merge so that ~ N^^^ electrons enter 
the island simultaneously. The origin of this effect is in 
the electron attraction mediated by the multi-polaronic 
effect associated with the confinement. A simple toy 
model of this effect was suggested in Ref. 



a) 



A 




FIG. 6. (a) Three point A, B, and C where the func- 
tion firo has to have extrema due to symmetry, (b) The 
oscillations of the thermodynamic potential for these three 
positions of the center. The minimum of these functions, 
shown by the bold line, represents the actual value of the 
total energy of the island. 

The dependence of the total energy on the number of 
particles can be obtained from our result using the small 
increment theorem, stating that the small corrections to 
all the thermodynamic potentials are the same:E]E3 



SE{N) = Sn [n{N)] . 



(3.40) 



This explains why in the dependence of the total energy 
on the number of particles we also observe a Kremlin wall 
structure similar to Fig. M. The amplitude and phase of 
the oscillations in Fig. ^ are in a very good agreement 
with those given by Eq. (3.37) minimized in the way 
shown in Fig. |^. 

We finally would like to notice that the appearance of 
a new terrace actually means starting a new crystalline 
row. This can be thought of as appearance of the pair 
of opposite dislocations on the boundary of the island. 
Thus the variations of the total energy discussed above 
are associated with the periodic formation of the defects 
in the island. This line of thinking will be further devel- 
oped for the case of the compressible crystal, where these 
defects are situated inside the island. 
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IV. COULOMB ISLAND: A NUMERICAL STUDY 

In this section we p resent the resuhs of numerical so- 
lution of p rob lem (1.2) with the unscreened Coulomb in- 
teraction (O). To obtain these results we employed the 
numerical technique described in the previous section. 
The total energy E{N) resulting from such a calculation 
was split into the smooth and the fluctuating components 
in the way described in the previous section. The smooth 
component was chosen to have the formO 



E (N) = (eV^)'/^ (7,iiv5/3+ 



(4.1) 



where r]i are some constants. The first term in this se- 
ries is the electrostatic energy. The next three terms are 
the correlation energy, the overscreening energy associ- 
ated with the screening of the external potential by the 
Wigner crystal, and the surface energy. The coefficients 
rji can be found from the best fit to the numerical data: 
= 6/5 • (37r/8)^/^ , r/2 = -1.0992 , r/g = -0.3520 , 
7]4 — 0.1499. The fluctuating part is displayed in Fig. |^a. 



a) 
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FIG. 7. (a) Fluctuating part of energy of the Wigner 
crystal island, (b) The number of particles adjacent to 
the center of the confinement. 

The curve in Fig. |^a has a quasi-periodic structure sim- 
ilar to the incompressible case (see Fig. ||). It consists of 
the sequence of interchanging deep and shallow minima. 
The positions of minima on this graph almost exactly 



coincide with those for the case of incompressible island. 
The amplitude of the oscillations does not change with 
N appreciably and is ~ O.le^/ Ka. The charging energy 
A(iV) experiences fluctuations « 15% correlated with the 
positions of the maxima and minima. 

The similarity of Figures |^ and |^a calls for the conclu- 
sion that the fluctuations in both cases are of the same 
origin. To test this idea we studied such a quantity as 
the number of electrons adjacent to the center of the con- 
finement. Those are defined as electrons nearest to the 
center, with dispersion in the distance to the center being 
less than 1 /3 of the lattice spacing. They represent the 
first crystalline shell closest to the center. The number 
of these particles is shown in Fig. 0b as a function of the 
total number of electrons in the island. The first thing 
to notice is the correlation of the latter graph with the 
fluctuations of energy in Fig. |^a. The deep minima in 
the energy are associated with having one electron next 
to the center. The shallow minima mostly correspond to 
having three electrons at the center. This in fact means 
that the center of the conflnement is situated in the cen- 
ter of the triangular crystalline face. Maxima of energy 
usually correspond to having two or four electrons next 
to the center. This is equivalent to saying that the cen- 
ter is in the middle of a crystalline bond. Thus we see 
that this correspondence is exactly the same as switch- 
ing between terms A, B, and C in the incompressible case 
(see Fig. |6|). This analogy will be further developed in 
Section | 

Let us now examine the conformations of electrons. 
Some of them are shown in Figures I and |. The first 
thing obvious from these Figures is that the surface of 
the island is not so rough as it was in the incompressible 
case (compare to Fig. At large N the surface contains 
no defects. The defects reside in the interior instead. 
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FIG. 8. The magic number configuration A*' = 85. Six 
five-coordinated particles associated with disclinations are 
marked by rings. Triangulation of the neighborhood of 
one of them is shown explicitly. 
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Let us define these defects. Elementary defect in a 2D 
lattice is disclination. The only possible form of such 
a defect in 2D crystals is the so-called wedge disclina- 
tion. It can be viewed as a wedge which is removed 
(inserted) from the crystal. The "charge" of the discli- 
nation is the angle, formed by the wedge. The mini- 
mum possible disclination charge for a triangular lattice 
is 7r/3. The disclinations can be identified with the par- 
ticles having an anomalous coordination number. In the 
triangular crystal the cores of the positive or negative 
disclinations are associated with the particles having five 
or seven nearest neighbors respectively, the normal co- 
ordination being six. Some examples of such defects are 
shown in Fig. |^. 

Dislocations are the pairs of positive and negative 
disclinations forming a dipole. They can be seen in Fig. ||. 
When number of electrons in the island is less than some 
critical value N « 150, there exist highly symmetric elec- 
tron configurations free of dislocations. These configu- 
rations can be realized only at some distinct values of 
N = N„i — 7,19,35,55,85, which we call magic num- 
bers. One of such magic number configurations N„i — 85 
is demonstrated in Fig. ||. On the other hand if iV > iV 
the dislocations are always present and the magic number 
configurations never exist. 




• VS/ • • • • ®V • 
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FIG. 9. The electron configuration for = 235. 
The particles having five and seven nearest neighbors are 
shown by the rings and the triangles respectively. 

We would like to stress again that all these defects are 
not observed directly on the surface. Moreover, if the 
number of electrons in the island is large, they do not 
come close to the surface. The same is true about the 
central region of the island. It is also usually free of de- 
fects. All the irregularities in the lattice are normally 
observed in a ring of width (2 ~ 3)a at a fixed distance 
from the edge. 

We would like to discuss now the accuracy of the re- 
sults we have obtained. We believe that the total energy 
for the cases shown in Fig. |^a is calculated with precision 
10^^(e^/K)^/'^A^^/'^. This was tested by reruns (starting 



from different initial conditions) and comparison with the 
results of simulated annealing. The configurations how- 
ever cannot be reproduced reliably. The thing is that 
several completely different configurations for the same 
number of electron can have very close energies, within 
10^^(e^/K)^/'^v4^^/'^ from each other. However the gen- 
eral features of all these low energy configurations de- 
scribed above remain true. Thus, although we cannot 
reproduce the position of every individual defect reliably, 
we believe that we can say that in general defects reside 
in a ring concentric to the surface of the width of a few 
lattice constants. 



We will see in Section VI that the interaction of defects 



is crucial in understanding of the energy fluctuations. In 
the next section wc discuss the general properties of the 
defect distributions. These properties are essential for 
the arguments given in Section VI. 



V. LATTICE DEFECTS IN THE CIRCULAR 
WIGNER CRYSTAL ISLAND 

When the triangular crystal is packed into the island 
of circular form the obvious incompatibility of these two 
structures leads to the appearance of the lattice defects: 
disclinations and dislocations. The important difference 
between the former and the latter is that disclinations are 
always present in the ground state. The number of discli- 
nations in the island is determined by the Euler's theorem 
(see below) and, hence, cannot be changed. Dislocations 
on the other hand may or may not appear depending on 
whether they are energetically favorable. Let us first ex- 
plain why the disclinations have to appear in the ground 
state. 

The disclinations within some region of lattice can be 
identified by the number of 7r/3 turns that one has to 
make to walk around it. Normally one has to make min- 
imum six such turns (consider for example the central 
point in Fig. ||). For the region having Nc such defects 
inside, the minimum number of 7r/3 turns is 6 — Nc. Look 
e.g. at the pentagon in Fig. It obviously contains one 
disclination, because the number of the aforementioned 
turns is five. Let us now mentally walk around the whole 
circular sample along its edge. As the surface of the is- 
land is circular we do not make any turns at all. Hence 
the total number of 7r/3 disclinations is exactly six: 



(5.1) 



In the simplest possible case shown in Fig. ^ there are 
only six disclinations in the sample. If both positive (re- 
moval of 7r/3 wedge) and negative (insertion of such a 
wedge) disclinations are present, the total disclination 
charge should be calculated accounting for the sign of 
these defects: 



Nr = N+ - n: 



(5.2) 



where and are the numbers of positive and neg- 
ative defects correspondingly. Such complex cases are 
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realized for example when there are a few dislocations in 
the island. Dislocation is a pair of positive and negative 
disclinations bound together to form a dipole. Hence ad- 
dition of a dislocation to the sample increases both iV+ 
and N~ keeping their difference the same. Consequently 
the number of dislocations is not con trol led by the topo- 
logical constraints expressed by Eq. (5.1). This equation 
can also be proved using Euler's theorem. The proof can 
be found in Appendix ^ 

Now let us turn to the dislocations. At zero tempera- 
ture they embody inelastic deformations. There are two 
main reasons for existence of such deformations: inho- 
mogeneity of the concentration of electrons and the pres- 
ence of disclinations. The first reason is not universal: 
the exact profile of density is determined by the form of 
confinement and interaction potential. The second one 
is universal, becaus e th e presence of disclinations is re- 
quired by theorem (5.1). 

The dislocations produced by the varying, crystal den- 
sity were recently discussed by Nazarov.E^I He noticed 
that the dislocation density must be equal to the gradi- 
ent of the reciprocal lattice constant: 



b{r) = a (r) z x Va ^ (r) . 



(5.3) 



Here b (r) is the density of Burgers vector, z is the 
unit vector perpendicular to the surface, and a = 



2/n(r)y/S is the varying in space lattice constant. 
This formula is easy to understand. First we notice that 
a dislocation adds an extra row to the lattice. Hence the 
density of dislocations is given by the rate of change of 
concentration of these planes equal to 1/a. To obtain the 
Burgers vector density we multiply this gradient by the 
elementary Burgers vector a (see also Ref. |lj) . 

This formula has an interesting implication for the 
question of applicability of the elasticity theory to the 
Wigner crystal. Consider a sample of size L. Assume 
that there is a variation of concentration in the sample 
induced by an external source Sn ~ UikU, where Uik is the 
strain tensor. Let us find the value of Uik at which dis- 
locations start to appear. This would imply that inelas- 
tic (plastic) deformations occur and the elasticity theory 
brakes down. The energy of purely elastic deformations 
is 



(5.4) 



Here Y is the Young's modulus. For inelastic deforma- 
tions we have 



ErN, 



(5.5) 



where Ec is the dislocation core energy and Nd is the 
total number of dislocations in the sample. Eq. ( |5.3| ) 
provides the following estimate for the latter: 



6n L 
—r- ^ u,k — - 
nLa a 



(5.6) 



Comparing the elastic and inelastic energies we conclude 
that the deformations are completely elastic when 



Uik 



< 



Ec 
YaL 



(5.7) 



For the Wigner crystal with the Coulomb interaction us- 
ing Ec ^ / Ka, Y ~ e^/^a^, and Uik ^ u/L, where u is 
the characteristic displacement, we obtain 

u<a. (5.8) 

Hence equilibrium elastic displacements cannot exceed 
the lattice spacing for the Wigner crystal. This agrees 
with the known result that the elasticity theory has a 
zero-Radius of convergence with respect to the strain ten- 
sor .Ej In effect we conclude that the radius of convergence 
depends on the spacial scale of the problem and in the 
mac roscopic limit L — > cxd is indeed zero, according to 
(^.7|) . This surprising result emphasizes the difference 
between electron crystal and crystal consisting of heavy 
particles (atoms). In an electron crystal the relaxation 
time provided by tunneling is smaller than the time of 
experiment and, hence, the system can find the ground 
state. In an atomic crystal the relaxation time is typ- 
ically much larger than the experimental time and the 
system can remain in the metastable state described by 
the elasticity theory. 

Below we compare the number of dislocations of two 
different origin. First we calculate the total number of 
dislocations produced by the non-uniformity of the elec- 
tron density. We consider the case of the Coulomb in- 
teraction. For this case the corresponding electrostatic 
problem (ignoring the discreteness of the charge of elec- 
trons) can be solved exactly.llj The solution for the den- 
sity of electrons can be shown to be a "hemisphere" : 



"W = "o\/l--^ 



( 



V 8Ak 



2\ 1/3 



(5.9) 



According to Eq. (5^) the density of dislocations induced 
by varying electron density is 



nd{r) 



h{r) 



V3i?2(i_r2/-R2) 



3/4 



(5.10) 



Here h = a\^/2 is the elementary Burgers vector, equal 
to the distance between crystalline rows. The total num- 
ber of dislocations is readily obtained by integrating this 
distribution: 



N, 



dl 



2TTrdrnd{r) w 4.08 • iV^/^ 



(5.11) 



The electrostatic formula (5.9) is correct if the number 
of electrons is large. In the opposite case of a small island 
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the density profile is very far from being jE hemisphere. 
In this case the density is almost constant .□ This is obvi- 
ous from looking at Figs. ^ and ^ above. To evaluate the 
variation of the electron density we will find the lattice 
constant on the surface of the island and compare it to 
that in the center. To this end we notice t hat n ear the 
edge the density of dislocations given by Eq. jslol ) would 
become infinitely large. However it cannot exceed the 
density of electron themselves, given by Eq. (5.9). The 
position of the edge of the crystal can therefore be deter- 
mined matching the dislocation density and the density 
of electrons. As a result we obtain for the lattice constant 
on the edge the following expression: 



(5.12) 



This expression agrees with conclusions of Ref. || obtained 
in a different way. Our numerical data agree very well 
with this theory and give the following coefhcient: 



A = 0.88ao/'i?i/5, 



(5.13) 



The interelectron distance on the edge given by this for- 
mula differs very slightly from the lattice constant in the 
center of the sample if the radius of the island is not too 
large. To evaluate the number of dislocations for a small 
island we approximate the electron density profile with a 
parabola: 



n{r) ~ no 



r 



1 _ m 
A2 



Applying Eq. (5.3) to this expression we obtain 



(5.14) 



(5.15) 



In the consideration below we will assume that the num- 
ber of electrons in the island is not too large. Therefore 
the electron density is almost uniform. 

As it was mentioned above, another reason for appear- 
ance of dislocations is the stress, produced by disclina- 
tions. The latter deform the lattice enormously. These 
deformations can be significantly reduced by introduc- 
ing dislocations into the lattice. This process is usually 
referred to as screening. The screening of disclinations 
by dislocations has been studied before in the context 
of the hexatic liquid - homogeneous liquid transition.Ea 
The screening is that case is accomplished by polar- 
ization of the thermally excited dislocations above the 
Kosterlitz-Thouless transition. It was therefore treated 
in the framework of the linear Debye-Huckel approxima- 
tion. Here we deal with the case when there are no ther- 
mally excited dislocations and all the Burgers vector den- 
sity necessary for screening exists due to the appearance 
of new dislocations. 

The phenomenon can be most easibL. understood con- 
sidering the total disclination dcnsity:lla 



where the first term is the density of free disclinations 
while the second is the density of disclinations induced 
by the varying in space density of dislocations (the lat- 
ter are the dipoles formed from the former). The elastic 
energy of the crystal can be written in terms of this total 
defect density: 



-4 



d^q |stot(g)r 

(27r)2 g4 



(5.17) 



For case of the 
Q,g2^3/2^^^ where a = 0.9804. 



where Y is thej— Young's modulus 
Coulomb crystaO Y 
We notice that for the large distances q ^ the numer- 
ator in this expression is pushed to zero by the q'^ term 
contained in the denominator. Thus, ignoring the effects 
at the distances r ~ a, one can write the condition of the 
perfect screening: 



Stot{r) = 0. 



(5.18) 



Stot(^") = s(r) - eiuVihkir), 



(5.16) 



It is instructive now to consider one disclination in the 
cente r of an infinite sample: s(r) — soS{r). Equations 
( ^.18 ) and ( 5.16| ) have in general many solutions. Two 
of them are easy to guess: bi^ = so/27rr, bir = and 
b2x = 0, b2y = SoQ{x)S(jj). They represent the Burg- 
ers vector rotating around the disclination and the solu- 
tion in the form of a grain boundary respectively. They 
are different by the longitudinal Burgers vector density: 
b± — b2 + V/, where / is some scalar function, and 
hence have the same elastic energy. This non-uniqueness 
i s a co nsequence of essentially mean- field character of Eq. 
( ^.18 ). Discreteness of the dislocations removes this de- 
generacy in favor of the grain boundary. The idea is 
that the fluctuations of the local elastic energy caused 
by the discreteness of the defects can be estimated as 
E ~ ya^ln((r) /a) per defect, where (r) is the average 
distance between them. In the former distribution bi 
this average distance grows with the size of the sample 
(r) ~ V La. In the case of the grain boundary the average 
distance is maintained of the order of the lattice spacing. 
Hence the self-energy of defects in the latter case does 
not grow with the size of the sample. 

The next logical step is to consider six disclinations 
in the circular crystalline sample. A possible solution 
for the distribution of defects is shown in Fig. |l^. It 
is prompted by the results of the numerical simulations 
showing that in the majority of cases the defects (disloca- 
tions and disclinations) are situated in a ring concentric 
to the surface of the island. The regions adjacent to the 
center and to the edge of the island appear to be free of 
defects. We assume that the disclinations form a figure 
close to a perfect hexagon. The dislocations form grain 
boundaries connecting the disclinations. This is done to 
smear out the charge of disclinations in accordance with 
the screening theory expressed by Eq. (5.18|). 
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FIG. 10. The distribution of defects in a compressible 
island. Rings and arrows show the positions of disclina- 
tions and the direction of the Burgers vector density, re- 
spectively. The gray disk shows the free of defects central 
region. 

Next we calculate the distance from the surface to the 
layer of defects. We consider the case of a uniform crys- 
tal. The central region of the island can be formed with 
no deformations in it. It can be thought of as a piece 
of the uniform crystal having a circular form, the same 
as the one considered in Sections || and [II . The central 
region is shown in Fig. |l^ by a gray disk. The region 
between the surface of the island and the layer of defects 
is free of defects too. It cannot be free of elastic deforma- 
tions however. These deformations are the same as those 
of a 2D rectangular elastic rod, two opposite edges of 
which are glued together. To find the optimum thickness 
of this rod w (see Fig. |l^) one has to balance the energy 
associated with this bending, which tends to reduce w, 
and the surface energy of the grain boundary. The lat- 
ter is proportional to the length of the grain boundary 
and therefore has a tendency to increase w. The bending 
energy can be calculated from the elasticity theorycEl 



12 R 



(5.19) 



The surface energy of the grain boundary can be esti- 
mated from the core energies of the dislocations: 



ErN, 



(5.20) 



The total number of dislocation follows from elementary 
geometric calculation done for the incompressible crys- 
talline region of radius r = R — w: 



Nd2 = 12(1 - Vi/2)r/h, 



(5.21) 



where h = a\^/2. Minimizing the total energy E ~ 
Eel + Eg\j we find the optimum at: 



w 



1.5 



EcR 
Ya 



(5.22) 



In the derivations of (5.22) we used Ec — O.lle^n^/^/K 



and Y = O.dSO^e^n^^yKiReL M). This resuh can be 



also obtained from our condition of stability of a crystal 
(^.8|) . Indeed the strain tensor in the ring can be esti- 
mated as ^ w/R. The corresponding displacement 
u ~ uP' / R cannot exceed the lattice spacing a, according 
to Eq. (5.8). The estimate for the width of the ring w ob - 



tained from this argument is consistent with Eq. ( 5.22 ). 

Now we would like to find the condition at which the 
dislocations stay on the surface of the island. We will 
consider the case of the short-range interaction ( |1.4D . To 
this end one has to compare the energy of the defects on 
the boundary with that in the bulk of the crystal. The 
former is associated with the roughness of the surface 
(see Fig. |). It can be estimated as a product of the total 
number of particles there a/ZV, the typical force acting on 
a particle on the surface F{R) = 2AR, and the character- 
istic deviation of the shape of the surface from the circle, 
given by the lattice constant: (5i?sur ^ • F{R) ■ a. 
The energy of the defects inside the island is given by 
(^^buik ^ ■ Ya^XriN. The first term in the right- 
hand side is the number of defects and the second one 
is the typical energy per dislocation. It is estimated as 
the typical interaction energy of two dislocations. Note 
that we neglect the core energy of dislocations as it is 
small compared to the Ya^ for the crystal with short- 
range interaction. Indeed the former is of the order of 
the correlation energy: ~ ?7o exp(— a/s). The Young's 
modulus for the considered system can be estimated as 
the second derivative of the interaction potential between 
two particles: 



Y ■ 



■ exp 



i-i) 



(5.23) 



It is clear then that Ec/Ya? ~ s^/a? <^ 1. Comparison 
of SEsu-r with (5i?buik gives the condition that the defects 
stay on the surface: 



(5.24) 



Y\nN > A^/N 



The coefficient A can be conveniently related to the inter- 
action strength by balancing the forces acting on a par- 
ticle on the surface: AR ^ Uo/s ■ ex p(— a /a). Combining 



this equation with Eq. ( 5.23 ) and ( 5.24 ) we obtain the 
following condition for the defects staying on the bound- 
ary: 



a In iV > s 



(5.25) 



This coincides with the condition of uniformity of the 
crystal (see Section HI). Note that this entire consid- 
eration is based on the assumption of uniformity of the 
crystal. 

Finally in this section we evaluate the critical number 
of electrons N* at which the number of dislocations due 
to the inhomogeneity of the electron density Ndi and 
produced by the screening of disclinations Nd2 b ecom e 
equal. To do so we compare Eq. ( 5.15 ) with E q. (5.21 ). 
They become equal if A « 1.55ao or, using Eq. ( ^.12 ) 



TV = AT* ~ 700. 



(5.26) 
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Therefore if the number of electrons in the island is 
smaller that N* we dislocations are mostly due to the 
screening of disclinations and are arranged into the grain 
boundary (see Fig. ^o|). In the opposite case the disloca- 
tions are generated in the interior according to Eq. (5.3). 

The above consideration is essentially mean-field, 
treating the density of dislocations as a continuous quan- 
tity. In the next section we give an argument that the 
discreteness of dislocations is responsible for the fluctua- 
tions of the elastic energy of the Wigner lattice. 



VI. ELASTIC BLOCKADE 

The number of dislocations in the Wigner crystal is- 
land is of t he or der of the number of crystalline rows in 
it [see Eq. (|^]: 



Nd - R/a ^ 



N. 



(6.1) 



Hence, it grows while the island in filled with electrons. 
On the average one dislocation appears as N increases by 
SN ~ y/N. The elastic energy stored due to the deviation 
of the number of dislocations from the average is relaxed 
when a new one is added. This phenomenon is simi- 
lar to the Coulomb blockade, with the words "electron" 
and "electrostatic energy" replaced by "dislocation" and 
"elastic energy" . 

Using this analogy to the Coulomb blockade one can 
easily estimate the order of the fluctuations of energy. In 
the Coulomb blockade case the fluctuations of the elec- 
trostatic energy are given by the following expression: 



SE : 



2C 



(6.2) 



where — 1/2 < SN < 1/2 is the deviation of the number 
of electrons in the dot from the average, determined by 
the gate voltage. Following the convention of the above 
mapping one has to replace e^/C by the characteristic 
energy of interaction of two dislocations in the island 
Uq ~ Ya^ ~ e^/na. One also has to replace SN by the 
deviation of the number of dislocations from the average 
SNd ^ SN/ \/N . As a result we obtain: 



Kfl N 



(6.3) 



Here a « 0.5 is a numerical constant. The approximate 
value of this constant is obtained from the numerical re- 



sults of Section As the total period of such an "elas- 
tic" blockade is Tat ~ VlV, the maximum fluctuation of 
energy that can be reached is (5i?max ~ ae^ jna. 

Another implication of this analogy is the appearance 
of the elastic blockade peaks. In the case of Coulomb 
bloc kade the intersection of two terms described by 
Eq. ( |6.2| ) for different numbers of electrons in the dot 
produce a spike in the conductance through the dot. 
At this point a new electron enters the dot. For the 



"elastic" blockade at the similar point a new disloca- 
tion enters the crystalline island. The chemical poten- 
tial of an electron has a discontinuity of the order of 



A/i ~ dEldN ae^/naT, 



N 



-ae /R. The charging 



energy at this point fluctuates by a value 
dS^ 



SA = 



dN 



1 



(6.4) 



which is of the order of the average A = fie^ / kR, fi be- 
ing another constant. This however does not lead to the 
bunching in the charging spectrum as it did in the case 
of the short-range interaction, because of the smallness 
of the numerical constant a. 

Let us consider the elastic blockade in more detail. In 
the previous section we derived the condition at which 
the crystal can be considered to be uniform. If the num- 
ber of electrons in the island is less than some critical 
value N* ~ 700, the variations in the density of electrons 
can be ignored and the number of dislocations associated 
with the variable density is small. Below we examine 
these two regimes separately. 

i) Almost uniform crystal: N < N* . 

As it follows from the previous section in this regime 
the density of electrons in the island is almost uniform. 
The crystal is packed into the circular form by forming 
two monocrystals: the almost circular internal region, 
free of elastic deformations, similar to the incompressible 
island; and the external ring (see Fig. |lo|). The inter- 
face between these two monocrystals is a grain boundary, 
which can be viewed as a string of dislocations. 

Assume now that the position of the center of the con- 
finement relative to the lattice ro is fixed. This coor- 
dinate has been introduced earlier in Section III. One 



can then calculate the energy of the island as a function 
of the number of electrons and the position of the cen- 
ter Er„(N). Due to the symmetry considerations given 



in Section III this function has extrema at the point A, 
B, and C shown in Fig. ^ similar to the incompressible 
island. Filling of the island results in switching among 
three energetic branches Ea{N), Eb{N), and Ec{N), 
every time choosing the lowest. Below in this subsection 
we estimate the correction to the total energy SEr^ (N) . 

The grain boundary has a tendency to be at a fixed 
distance from the center r — R — w, calculated in Sec- 
tion If position of the center relative to the crystal is 
fixed, filling of the island brings about expanding of the 
mean grain boundary position with respect to the crystal. 
Hence, periodically the grain boundary has to intersect a 
new crystalline row. At this moment in the correspond- 
ing incompressible problem a new terrace appears on the 
boundary of the island. In the Coulomb problem the 
grain boundary is submerged into the island. Therefore 
at this moment in the compressible case a new row is 
added to the island. Since the termination points of the 
new row can be thought of as dislocations, we can say 
that a new pair of opposite dislocations appears on the 
grain boundary (see Fig. with the gray internal region 
shown in more detail in Fig. |^). 
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The corresponding fluctuation of energy can be eval- 
uated as the energy of interaction of two opposite dislo- 
cations having charge SNj, = Nd — Nd, where Nd is the 
actual number of dislocations in the neighborhood of the 
ne"W-.row, Nd is the average one. This energy is given 



(6.5) 



The typical distance between defects / is given by 
Eq. $A^ : I ~ VR-a. The average number of dis- 
locations is given by the number of crystalline rows 
Nd = '^ii'f — £ifo)/h^ e-i being the unit vector normal 
to the series of rows considered. Therefore 



6Nd = J2 



r - eiVo 



(6.6) 



where by {• • •} we assume taking the fractional part. The 
oscill ation s of energy are given by the formula similar to 
Eq. (^!37|) 



SEro = ^ SEo{r ~ e,ro), 



where 

SEoir) 



{1} 



1 

12 



IniV 



(6.7) 



(6.8) 



are the fluctuations of energy given by (6.5). The last 
term in the square brackets is chosen so that SEQ{r) = 0. 

The oscillations of energy calculated in this way for 
three positions of the center Tq = A, B, and C are shown 
in Fig. |ll|. 
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N 



FIG. 11. The fluctuating part of the energy of the 
Wigner crystal island calculated from the interaction of 
dislocations. 

As usually the lowest one is to be chosen. Compari- 
son to the numerical experiment (Fig. |^) shows, that this 
model predicts very well the phase and the shape of the 
oscillations. The amplitude however is smaller by a factor 



of 2 2.5 than observed in the simulations. This differ- 
ence can be explained by the influence of the boundary 
of th e sample on the interaction energy of dislocations 
( |6.5|) . Indeed the boundary of the island is at the dis- 
tance w « ■ a from the ring of dislocations. This 
scale is smaller than th e char acteristic distance between 
dislocations I « {2^3)VR ■ a. Hence the boundary of the 
island should play an important role in the interaction of 
dislocations. 

To understand the character of the correction we no- 
tice that the boundary in the numerical experiment is 
almost not deformed. The island tends to preserve its 
circular shape. Hence we can assume that the normal 
displacement of the boundary is zero. Assume now that 
the dislocations are much farther away from each other 
than from the boundary I ^ w. In practice we have 
I ~ 5w. The problem of interaction of the pair of dis- 
locations in the presence of the boundary can then be 
solved using the method of images. It is easy to see 
that at the zero normal displacement boundary condi- 
tion a dislocation with the Burgers vector perpendicular 
to the boundary produces an image of the same sign. In- 
deed such a dislocation is a crystalline semi-row, parallel 
to the surface. The zero normal displacement boundary 
condition can be realized by putting a parallel semi-row 
on the other side of the boundary. Hence effectively the 
boundary doubles the charge of each dislocation. The 
interaction energy, being proportional to the square of 
charge, would acquire a factor of four due to such im- 
ages. However as the crystal exists only in the semi-space 
the deformations only there contribute to the elastic en- 
ergy. Therefore the overall factor that arises due to the 
presence of the boundary is 4 • 1/2 = 2: 



Eu 



ndary(f) = 2Efrcc{r) 



2tt 



In 



id 



(6.9) 



This factor can explain the discrepancy between our cal- 
culation and the numerical experiment. 

ii) Non-uniform crystal: N ^ N* . 

In this case as it follows from Eq. ( |3.10 ) there are 
Nd ~ -v/iV dislocations in the island associated wit h the 
variations of the electron density. According to Eq. ( 5.10] ) 
these defects have to be present in the bulk of the crystal. 
Hence they cannot be arranged into grain boundaries as 
in the case considered above. The nearest dislocations 
repel each other and, one could imagine that they form a 
crystal themselves.t2l We argue however that this cannot 
be the case. The reason is that dislocation can occupy 
a fixed position within an electron lattice cell. Hence 
the incommensurability of the electron and dislocation 
lattices eventually produces frustration and destroys the 
long-range order in the dislocation lattice. Therefore we 
expect that the dislocations form a glassy state. 

This implies that the coherence of the crystalline rows 
in the electron crystal itself is destroyed. We think there- 
fore that in this regime the variations of the total energy 
are not periodic. They retain however the general fea- 
tures described in the introduction to this section. To 
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understand how these features arise in this particular 
case we would like to study the reconstructions in the 
dislocation lattice caused by addition of a new electron. 
To this end it is necessary to consider two cases: when 
this addition docs not change the total number of de- 
fects, and when the number of defects is increased. The 
former case takes place most of the time, while the latter 
happens once in \/7V electrons when a new defect has to 
appear. 

Consider the first case. Electron can be added into 
the core of dislocation increasing the length of the extra 
crystalline row associated with it. The dislocation after 
this is displaced by the lattice constant. Let us find the 
maximum change in the elastic energy of the island as- 
sociated with such a displacement. It can be expressed 
through the interaction between dislocations: 



(PU{r) 



(6.10) 



where U{r) = Yh? ln(r)/47r is the dislocation interaction 
energy. The reason why the second derivative is relevant 
to the calculation of this energy is that the dislocation is 
in equilibrium before adding the new electron. Hence the 
gradient of the self-energy of the dislocation is zero. Now 
we have to remember that there are ^/N dislocations in 
the island. Hence the actual addition energies range from 
zero to ^£^max with the average level spacing 



5E 



6E^ 



N 



e a 
1^' 



(6.11) 



This quantity describes the energy one pays when adding 
an electron at the fixed number of dislocations. This is a 
correction to the energy spacing. Multiplied by the total 
number of electrons between two consecutive additions 
of dislocations ^/N it gives the variation of the chemi- 



cal potential 5E^/N ~ / kR. This variation should be 
equal to the discontinuity of the chemical potential when 
a new dislocation is added to the island, as the average 
correction to this quantity due to elastic effects does not 
grow. Thus the drop of chemical potent ial is equal to 
~ / kR. This is consistent with Eq. (^.4|). 
Thus elongated crystalline lines draw dislocations from 
the center to the periphery of the island. Eventually a 
new row has to be inserted in the center. This event 
manifests itself in the appearance of a new pair of dis- 
locations. The distance from the center ^ at which they 
are situated can be found from Eq. (5.1C) by stating that 



This results in 



(6.12) 



(6.13) 



The energetics of switching between two branches corre- 
sponding to different number o f di slocations has already 
been discussed above [see Eq. (|6.4[)]. 



VII. CONCLUSIONS 



In this section we would like to discuss the degree of 
universality of our results. The first question is what 
happens if confinement potential V{r) is not parabolic. 
The answer is almost obvious for the hard-disc interac- 
tion. In this case electrons are added to the crystal on 
the equipotential determined by the level of chemical po- 
tential. Therefore only the shape of this equipotential 
and the confinement potential gradient matter for en- 
ergy fluctuations. It is obvious that situation is almost 
identical for any isotropic confinement V{r). It is also 
easy to generalize our calculations for the confinements 
of oval shapes. 



In the latter case the theory developed in Section ^1 
can be applied with a few modifications. As it was ex- 
plained the variations of energy are associated with pe- 
riodic intersections of the equipotentials with the crys- 
talline rows. The main contribution to these oscilla tions 
comes from the "coherence" spots of the size ~ ^ R - a 
tangential to the rows. The oscillations of energy are 
sensitive then only to the properties of the equipoten- 
tial in the immediate vicinity to those spots. If it can 
be well approximated there by a circle Eq. (3.37) should 
hold with additional phase shifts introduced into the ar- 
guments of the contributions from different rows. These 
phase shifts arise due to the deviation of the global shape 
of the equipotential from the circle, and express the inco- 
herence of contributions coming from six different series 
of rows. In addition to that each contribution should 
be rescaled according to (x y/RlT^^ where Ri is the 
curvature radius of the corresponding equipotential. We 
would like to emphasize again that this theory works only 
for an island of oval shape. It is not applicable for ex- 
ample to a square. All of the curvature radii Rt have 
to be of the order of the size of the island. In general 
as the contributions from different terraces are not co- 
herent anymore, the overall amplitude of oscillations has 
to be smaller in the oval case compared to the case of 
rotational symmetry. 

Let us discuss the universality of the results with re- 
spect to the choice of the interaction potential. As we 
have seen both short-range and Coulomb interactions 
give rise to the fluctuations of energy of the same func- 
tional shape. We argue that the other forms of long-range 
interactions, logarithmic for instance, bring about simi- 
lar results. Thus the energy of a disorder-free cylinder 
or disk of a type-II superconductor filled with the fluxoid 
lattice or a rotating cylindrical vessel of the superfluid 
helium as a function of the number of vorticesEil should 
experience oscillations similar to those in Fig. 0. To check 
this prediction we have performed a numerical calcula- 
tion in our model with the interaction U(r) = Uq ln(l/r). 
It revealed a quasi-periodic correction to the energy of the 
same functional shape as shown in Fig. ^ The amplitude 
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of the correction was consistent with the conclusions of 
Section Vl. The configurations were also very similar 
to the observed in the Coulomb case. In particular the 
correlation between the position of the center of the con- 
finement relative to the crystal and the oscillations of the 
energy is also observed. 

In conclusion we have studied the charging spectrum 
of the crystal formed by particles in the parabolic con- 
finement. We considered two forms of the interactions 
between particles: the short-range and Coulomb interac- 
tions. In the computer simulations employing the genetic 
algorithm we have observed the oscillations of the ground 
state energy which have a universal form, independent of 
the form of interaction. We attribute these oscillations to 
the combination of two effects: periodic additions of new 
crystalline rows and hopping of the center of the confine- 
ment relative to the crystal. The hops are separated by 
addition of ~ N'^^'^ electrons. These hops make a dra- 
matic difference for the addition spectrum of the island. 
In the case of the short-range interaction they make the 
charging energy negative, so that ~ N^/'^ new electrons 
enter the island simultaneously. This apparent attraction 
between electrons is a result of the confinement polaron 
effect discussed in Ref. In the case of Coulomb inter- 
action such a hop results in an abrupt « 15% decrease 
in the charging energy. 



where v, /, and e are the numbers of vertices, faces, and 
edges contained in the graph formed by our triangula- 
tion. All the vertices of the graph can be divided into 
the internal ones Vi, belonging to the bulk, and the ones 
on the surface Ve- The same can be done for the edges: 



C — Gi ~\- Cg . 



(A2) 



As all the faces of our figure are triangular by construc- 
tion the following relationship, expressing the general 
balance of edges, is true: 



3/ = Ge -1- 2ei. 



(A3) 



Next we can relate the number of edges to the number 
of vertices. As it was mentioned above the bulk vertices 
are connected to six edges while the surface ones to four. 
The exceptions are the cores of disclinations. They have 
an anomalous coordination number. Expressing now the 
overall balance of edges one can write: 



2e — 6vi — Svi + — 6ve 



(A4) 



where Svi and Sve are total deviations from the normal 
coordination numbers for the internal and external ver- 
tices respectively. One can finally use the obvious fact 
that 
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APPENDIX A: THE TOTAL DISCLINATION 
CHARGE OF THE ISLAND AND THE EULER'S 
THEOREM 



(A5) 



to solve the system of equations 
tain: 



Nc = Svi + SVe = 6. 



W - (A4) and to ob- 



(A6) 



This completes the proof of Eq. (5.1). We would like to 

kill 



notice that a similar theorem is wel 
gular crystal on the surface of a sphere 

Nr = 12. 



for a trian- 



(A7) 



To prove Eq. ( p.lD one has to first consider some tri- 
angulation of the electron lattice. It is convenient to 
consider the triangulation in which every electron is con- 
nected by edges to the nearest neighbors. In the case of 
triangular lattice an electron in the bulk has six nearest 
neighbors, while an electron on the surface of the sample 
has only four. For some electrons however this number 
can be different. For example, an electron in the core 
of 7r/3 disclination has only five nearest neighbors (see 
Fig. H) . Electrons on the surface can have the coordina- 
tion number equal to three. We associate such electrons 
with the 7r/3 disclinations stuck to the surface. 

Let us now use the Euler's theorem. For our case it 
says: 
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